Zone-Center Dynamical Matrix in Magnetoelectrics 
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In ordinary dielectrics the dynamical matrix at the zone center in general is a nonanalytic function 
of degree zero in the wavevector q. Its expression (for a crystal of arbitrary symmetry) is well 
known and is routinely implemented in first principle calculations. The nonanalytic behavior occurs 
in polar crystals and owes to the coupling of the macroscopic electric field E to the lattice. In 
magnetoelectric crystals both electric and magnetic fields, E and H, are coupled to the lattice, 
formally on equal footing. We provide the general expression for the zone center dynamical matrix 
in a magnetoelectric, where the E and H couplings are accounted for in a symmetric way. As in the 
ordinary case, the dynamical matrix is a nonanalytic function of degree zero in q, and is exact in the 
harmonic approximation. For the sake of completeness, we address other issues, and in particular 
we solve a problem which might arise in first-principle implementations, where — differently than 
here — the basic fields are E and B (not H). 

PACS numbers: 63.20.-e, 75.85.+t, 63.20.dk 
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I. INTRODUCTION 

The zone-center dynamical matrix in polar dielectrics 
is comprised of an analytic term and a nonanalytic term; 
the latter accounts for the long range of Coulomb inter- 
actions, or equivalently for the coupling of macroscopic 
electric fields with the lattice. In high symmetry situa- 
tions the nonanalytic term is responsible for the famil- 
iar longitudinal-transverse splitting of zone-center opti- 
cal modes. The explicit form for the zone-center dynam- 
ical matrix in crystals of arbitrary symmetry was first 
provided in 1962 by Cochran and Cowley. 1 Their phe- 
nomenological formula is exact within the harmonic ap- 
proximation and has been implemented much later in 
some first-principle codes. 2-4 

Magnetoelectrics (MEs) are insulators where electric 
fields control magnetization, and conversely magnetic 
fields control polarization; they attracted considerable 
theoretical and technological interest in recent times. 5-16 
In such materials both fields (electric and magnetic) 
are coupled to the lattice: therefore both contribute 
to the nonanalytic term and (in high symmetry cases) 
to the longitudinal-transverse splitting. In Ref. 15 it 
is shown that the time-honored Lyddane-Sachs- Teller 
relationship, 17-20 which applies to ordinary dielectrics 
only, generalizes in a perspicuous way to MEs. In the 
present paper we provide the exact form of the zone- 
center dynamical matrix in a crystal of arbitrary symme- 
try, where the nonanalytic term accounts for the coupling 
of both fields to the lattice in a very symmetric way. 

So far, we have not specified which macroscopic fields. 
In lattice dynamics, the natural choice is the pair (E, H). 
In fact, they are both longitudinal, while the fields D and 
B are both transverse. 21 ' 22 Because of this key feature the 
contribution of (E, H) to the restoring forces appears as 
a nonanalytic term, while the analytic one corresponds to 
setting E = H = 0. In contrast to this situation, the nat- 
ural choice in the framework of first-principle calculations 
is the pair (E, B); in fact a second order expansion of the 



energy per cell with the ordinary periodic boundary con- 
ditions correspond to E = B = and, therefore, does not 
provide as such the analytic term in the force-constant 
matrix. This problem is addressed and its solution given. 

Another issue addressed in the present work is the mi- 
croscopic nature of the coupling of magnetic fields to the 
lattice, which at first may appear as counterintuitive. In 
fact, a microscopic magnetic field does not exert any force 
on a nucleus at rest. 

The plan of the paper is as follows. In Sec II we estab- 
lish our notations and we review the Cochran-Cowley for- 
mula for ordinary dielectrics. 1,2 in Sec. Ill we show how 
this formula generalizes to MEs, arriving at our major re- 
sult, Eq. (25). We then deal with other related issues. In 
Sec IV we discuss the relationship of the present work to 
a recent generalization (to MEs) of the Lyddane- Sachs- 
Teller relationship. 15 In Sec V we show how a magnetic 
field in MEs does indeed exert a force on a nucleus at 
rest. In Sec VI we address the H versus B issue, cru- 
cial for any first-principle implementation of the present 
result. Finally, in Sec. VII we draw our conclusions. 



II. ORDINARY DIELECTRICS 



For the sake of completeness, as well as for estab- 
lishing our notations, we provide here a derivation of 
the Cochran-Cowley 1 phenomenological formula for the 
zone-center dynamical matrix in ordinary dielectrics of 
arbitrary symmetry. Whenever the crystal is polar and 
insulating, the dynamical matrix shows a nonanalytic be- 
havior at the zone center, which accounts for field-lattice 
coupling. Only electrical fields are considered in this Sec- 
tion. 
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A. Notations 

We use compact notations leaving the Cartesian in- 
dices implicit throughout. The electronic (called also 
clamped- nuclei or "static high frequency") dielectric 
tensor 23 is real symmetric and indicated &s Sqo \ the zone- 
center analytic part of the force-constant matrix, indi- 
cated as C ss i is real symmetric for a simultaneous ex- 
change of both its (implicit) Cartesian indices and its ba- 
sis indices s,s'. This analytic part yields by definition 
the second order expansion of the energy in the lattice- 
periodical displacements u s at E = 0. Whenever the 
force-constant matrix is computed from first principles, 
the ordinary choice of periodic boundary conditions is 
equivalent to assume E = 0; the computation provides 
then C ss > directly. The magnetic analogue is different 
in this respect: a discussion is provided in Sec. VI. The 
Born-charge tensors Z* and Z* 1 * (transpose) are nonsym- 
metrical Cartesian tensors. We further define the unit 
vectors in the q-direction as 



q = 




q 1 = ( Qx % q z ) 



(1) 



Therefore the norm is q' q = 1, while the dyadic product 
qqf is the projector in the direction of q. 



B. D and E fields 

In presence of a long wavelength phonon of wavevector 
q, the solid is macroscopically homogeneous in the plane 
normal to q, while all macroscopic properties display a 
modulation in the direction of q. It is immediate to re- 
alize that the component of D(q) parallel to q and the 
component of E(q) normal to q both vanish: 21,22 



qtD(q) = 0, (l-qqt)E(q) = 0. 



(2) 



Whenever nonvanishing, both D(q) and E(q) are nonan- 
alytic functions of order zero in q. In the following equa- 
tions, it is tacitly assumed that only the leading term in 
q is considered. 

If P(q) is the macroscopic polarization due to a phonon 
in zero E field (transverse polarization), the fields D and 
E are related by 



D(q) = eoo E(q)+47rP(q). 
We now exploit Eq. (2) as follows: 

= q f D(q) = q t £ co E(q) + 4 7 rq t P(q) 

E(q)=qq+E(q). 
From these it easily follows that 

= q t £ 0O qq t E(q)+47rq t P(q) 



E(q) = 



47T 



q^ooq 



qq f P(q), 



(3) 



(4) 
(5) 

(6) 
(7) 



which can be interpreted as the depolarization field for 
an arbitrary q-direction. 



C. Equations of motion 

The equations of motion in the harmonic approxima- 
tion are 

-MV(q)u s (q) = / s (q), (8) 
where the forces and the fields, to leading order in q, are 

/ s (q) = -^C ss ,u s ,(q) + ^; t E(q) 

s' 

4-7T 

D(q) = -]>>> s (q)+ eoc E(q). (9) 

s 

In Eqs. (33) and (9) u s are sublattice displacements, 
M s the corresponding nuclear masses, and is the cell 
volume. From the second line of Eq. (9) it is clear 
that the phonon polarization in zero E field is P(q) = 
Z*u s (q), and the corresponding depolarization 
field at arbitray q is then 



E(q) = 



47T 



£1 qtfTooq 



q^qtz;u s (q) 



(10) 



The first line of Eq. (9) yields therefore 

4tt (Z*tq) (qtZ,*,) 



/ S (q) = -E 



c s 



q f £ooq 



Mq). (11) 



The quantitiy in parenthesis is indeed the usual expres- 
sion for the force-constant matrix at the zone center, in- 
cluding the nonanalytic term, first obtained in 1962 by 
Cochran and Cowley 1 and implemented much later in 
some first-principle codes. 2-4 This confirms that the ma- 
trix elements at the zone center are indeed nonanalytic 
functions, homogeneous of degree zero in q; we also re- 
mind that Eq. (11) applies to crystals of any symmetry. 

It is also expedient to provide an expression equivalent 
to Eq. (11) by restoring in it the normalized wavectors 
and indicating with 'P(q) = qq^ the projector in the q- 
dircction. Then 



/ s (q) = -£ 



C ss i + 



4tt z;tp(q)z;, 



£1 qteooq 



IV (q). (12) 



It is immediate to verify in either Eq. (11) or Eq. (12) 
that even the nonanalytic term is symmetric for a simul- 
taneous exchange of both the Cartesian indices and the 
basis indices. 

The nonanalytic term in the force constants accounts 
for an additional restoring force due to the depolariz- 
ing field; when the corresponding dynamical matrix is 
evaluated, this term can be interpreted as a generalized 
(squared) plasma frequency. We also notice that, at vari- 
ance with the familiar high symmetry cases, in a low sym- 
metry crystal the two terms (analytic and nonanalytic) 
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in the dynamical matrix do not need to commute. When 
this is the case, all zone-center modes are coupled to a 
nonvanishing E field (i.e. they are infrared active) while 
the C ss i by themselves do not correspond to any physical 
mode. 



III. MAGNETOELECTRICS 

In any linear ME the role of the 3x3 Cartesian tensor 
£oo is played by a 6 x 6 response matrix — called 1Z here — 
which yields the macroscopic fields (D, B) in terms of 
(E, H). We define the purely electronic (clamped- nuclei) 
response as 15 




where fi^ , and are the clamped-nuclei magnetic per- 
meability and ME coupling tensor, respectively. 

A. Depolarization and demagnetization fields 

If P(q) and M(q) are the macroscopic polarization 
and magnetization due to a phonon at E = H = 0, the 
analogue of Eq. (3) is 




We then exploit Eq. (2) together with its magnetic ana- 
logue 

qtB(q)=0, (1 - qqt)H(q) = 0, (15) 

to obtain the two scalar equations 

= q t e oqq t E(q)+qt aoc qqtH(q)+47rqtp(q)(16) 
= q t at c qqtE(q)+qt Aloc qqtH(q)+47rqtM(q). 

The scalars q^E(q) and qtH(q) are easily obtained by 
inverting the 2x2 matrix 

*<«-(&eS&J) (17) 

Finally we exploit once more E(q) = qq^E(q), H(q) = 
qq^H(q), to obtain both depolarization and demagneti- 
zation fields as 

(§$)—*-'<«(«)■ ci.) 

Both fields are longitudinal (parallel to q); Eq. (19) 
clearly generalizes Eq. (7) to the ME case, and coincides 



with it for a magnetically inert material. A comment 
about notations: in order not to overburden notations, 
we indicate with .M _1 (q) in Eq. (19) the 6x6 matrix, di- 
agonal over its Cartesian indices, whose three 2x2 blocks 
are indeed the inverse of M(q) in Eq. (17). 

Not surprisingly, the fields E(q) and H(q) for a given 
q direction depend on both P(q) and M(q), defined as 
the macroscopic polarization and magnetization in zero 
fields. This, indeed, is the hallmark of ME materials 

B. Equations of motion 

One starts from the most general expression for the 
free energy .F({u s }, E, H), where the sublattice displace- 
ments and the fields E, H are chosen as independent 
variables. ' 15 The coefficients of its second order expan- 
sion, entering the equation of motion, are tensorial mate- 
rial constants which in general are all independent from 
each other. The second derivatives at zero displacements 
are proportional to IZoc, while the mixed second deriva- 
tives of F with respect to the displacements and either 
E or H are the lattice coupling tensors (electrical and 
magnetic, respectively). 

In the present notations the forces and the fields, to 
leading order in q, are 

/ s (q) = -J2 CWMq) + ZfE(q) + C f H(q) (20) 




where £*, first introduced by Ihiguez, 7 is the magnetic 
analogue of the Born effective-charge tensor. Notice that 
in Eq. (20) the analytical term in the force-constants 
matrix C ss ' coincides by definition with the full force- 
constants matrix in zero E and H fields (not B). This 
feature deserves discussion, provided in Sec. VI. 

The macroscopic fields associated to a long- wavelength 
phonon of wavevector q and lattice displacements u s (q) 
is therefore, according to Eq. (19), 

( ) ~s < 21 » 

We simplify a bit the notations, by indicating from now 
on with "P(q) the "double projector", i.e. the block- 
diagonal matrix 

"W-^'w)- (22) 

We notice that 'P(q) is a 6x6 matrix, diagonal on the field 
indices, while we remind that A / J _1 (q), also a 6 x 6 ma- 
trix, is instead diagonal on the Cartesian indices. There- 
fore the two matrices 'P(q) and A / I _1 (q) commute. It is 
also expedient to define the ME lattice-coupling matrix: 

Z$ = (Z$,{t). (23) 
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Notice that Z\ is a 3 x 6 matrix (and Z s is 6 x 3). 
In these compact notations Eq. (21) becomes 



E(q) 
H(q) 



4-7T 

■- M-\ci)V(ci)J2z>s(<lh (24) 



and hnally replacing into the first line of Eq. (20) we get 
the generalized form of Eqs. (11) and (12) as 



/.(q) = -£ 



c s 



47T 

77 



'(q)- 

(25) 

Owing to the fact that 'P(q) and M 1 (q) commute, both 
terms in Eq. (25) are — as they must be — symmetric for 
a simultaneous exchange of both the Cartesian indices 
and the basis indices. In a low-symmetry crystal, the 
two terms do not in general commute, as indeed in the 
purely electrical case. Eq. (25) is the central result of 
this work; it is a consequence of the fact that, in ME 
crystals, both macroscopic fields E and H are coupled 
to long-wavelength modes, and both are therefore at the 
root of the nonanalytic term. 



IV. LYDDANE-SACHS-TELLER 
RELATIONSHIP 

We address in this Section only crystals whose symme- 
try is orthorombic or higher, in which case all crystalline 
tensors can be simultaneously diagonalized. This is e.g. 
the case for the paradigmatic ME crystal Cr203. If we 
choose q along a principal axis, then the zone-center op- 
tical modes are either longitudinal (u s parallel to q), or 
transverse (u s normal to q). 

Let us consider first, both for simple dielectrics and for 
MEs, the analytic term only in the zone-center dynamical 
matrix 



T-j(analytic) _ 



1 



--C s 



(26) 



When the Cartesian axes coincide with the principal 
axes, this matrix factorizes; for a crystal with an (N+ 1)- 
atom basis, each of the blocks leads to one zero eigenvalue 
(acoustic mode) and N nonzero eigenvalues, which we 
call w 2 and correspond to the (squared) optic frequen- 
cies of the system for a long wavelength-mode with q 
normal to the chosen principal axis. The restoring forces 
responsible for these modes do not have any contribution 
from the macroscopic fields (either E or H). 

If we consider instead long- wavelength modes whose q 
vector is parallel to the chosen axis, then the projector 
V(q) in both Eq. (12) or Eq. (25) act like the identity 
(either 1 x 1 or 2 x 2, respectively), and the dynamical 
matrix has (in general) different eigenvalues and eigen- 
vectors from the previous case. We indicate these longi- 
tudinal eigenvalues as cj 2 . 

In this work we have not yet addressed the response of 
the system to a genuinely static perturbation. In macro- 
scopic fields (either E or H) the total polarization and 



magnetization are due to the response of both the elec- 
tronic system and the lattice. 

For ordinary dielectrics, is it customary to indi- 
cate as £oo the electronic dielectric tensor and with 
£o the genuinely static one, which includes the lattice 
contribution. 23 In an high-symmetry situation we may 
choose a principal axis and consider the scalar eo/soo > 1 • 
The Lyddane-Sachs- Teller (LST) relationship in its orig- 
inal form 17 applies to a binary crystal, where the number 
of optic modes is N = 1; it states that eo/£oo = w 2 /w 2 , 
where the modes are polarized along the same chosen 
axis. Whenever N > 1, the ratio eo/^oo is related to 
a function of the w 2 and cj 2 . One of the expressions 
for this function has the form of the ratio of a weighted 
harmonic average of the a) 2 , over a weighted harmonic 
average of the w 2 . Obviously this yields the original LST 
relationship for N = 1. 

The LST relationship — either its original or general- 
ized form 1 ' 18-20 — is exact in the harmonic approxima- 
tion. The two members of the identity are greater than 
one in polar crystals, as a consequence of the same phys- 
ical effect: the coupling of the lattice with macroscopic 
electric fields. 

In MEs we may generalize Eq. (13) by including the 
lattice contribution 



Ho 



e a 
a Mo 



(27) 



When projected over a principal axis, both Tvloo and TZq 
become 2x2 matrices. It has been recently shown 15 that 
a scalar function of these two matrices takes the role of 
£o /^oo in a generalized LST relationship and equates the 
ratio of a weighted harmonic average of the w 2 , over a 
weighted harmonic average of the w 2 . At the root of the 
need for a generalization of the LST relationship is the 
fact that both fields E and H are coupled to the lattice 
on equal footing. The explicit form of the generalized 
LST relationship for ME crystals is given in Ref. 15, 
where the noncrystalline and/or anharmonic cases are 
also addressed. 



V. MICROSCOPIC ORIGIN OF THE 
FIELD-LATTICE COUPLINGS 

By definition, the Born-effective charge tensor 
yields the force exerted on the s-th nucleus at equilibrium 
by a macroscopic field E. This clearly appears from the 
first line of Eq. (9) or even Eq. (20) with H = 0. Equiv- 
alcntly its transpose yields the macroscopic polarization 
induced by a sublattice displacement u s at zero E field. 
According to the second line of Eq. (9) such polarization 
is in fact given by Z*u s /£l. 

The coupling tensor plays a similar role in the mag- 
netic case. If we set E = in Eq. (20), in its first line 
yields the force exerted on the s-th nucleus at equilib- 
rium by a macroscopic field H, while (* in its second line 



■5 



yields the macroscopic magnetization £*u s /f2 induced by 
a sublattice displacement u s at E = H = 0. 

The dual view just presented owes to the fact that both 
tensors are the second mixed derivatives, with respect to 
cither E or H and to sublattice displacements u s , of a 
free energy, where u s , E, and H are the independent 
variables. 7 ' 15 

While we have addressed macroscopic fields so far, we 
switch to microscopic fields next. The force f s on a nu- 
cleus of (bare) charge Z s at zero displacement is equal to 
Z S ~E S , where E s is the microscopic field at site s. Such 
statement is no longer correct in a pseudopotential frame- 
work; to keep things simple, we adopt an all-electron view 
throughout this Section. We remind that, by definition, 
the macroscopic field in the lattice-periodical case is the 
cell average of the microscopic one. 

In ordinary dielectrics we have therefore Zf~E = Z S E S , 
or equivalently 

E s - ^ZfE, (28) 

meaning that the tensor Zp/Z s yields the microscopic 
field at site s as a linear function of the macroscopic 
field. Notice that, while Z s is always a positive integer, 
the Born tensors fulfill the acoustic sum rule, and their 
sum vanishes: in the diagonal case they are real numbers 
bearing either sign. 

In MEs the force f s on the sth-nucleus at equilibrium 
is in general nonzero even at E = 0, provided that H ^ 
0, as shown by the first line of Eq. (20). Given that a 
magnetic field does not exert any force on a charge at 
rest, one may wonder about the microscopic origin of 
this force and of the corresponding coupling tensors. 
The explanation is that the force on a nucleus is always 
given by f s = Z S E S , but the microscopic field E s at the 
nuclear site is a linear function of both macroscopic fields 
E and H. In particular, E s is in general nonzero (and 
linear in H) even when E = 0. 

VI. FIRST-PRINCIPLE CALCULATION OF 
THE FORCE-CONSTANT MATRIX: H VERSUS B 

All of the above results — and in particular the expres- 
sions for the zone-center force constants, Eqs. (II) and 
(25) — are exact in the harmonic approximation, and ap- 
ply to either empirical models or first-principle calcu- 
lations. Indeed, the first occurrence of Eq. (11), due 
to Cochran and Cowley, 1 predates first-principle calcu- 
lations by almost three decades. 

For many dielectric materials the dynamical matrix 
is nowadays routinely computed from first principles, 2 
and for instance is provided by some codes in the public 
domain. 3 ' 4 At the zone center, the force-constant ma- 
trix has the form given by Eqs. (11) and (12). Therein, 
the ingredients of the second (nonanalytic) term are the 
clamped-nuclei dielectric tensor e m and the Born effec- 
tive charge tensors Z*: both quantities are usually com- 



puted using a linear response algorithm. As for the an- 
alytic part of the force-constant matrix C ss i , it is com- 
puted via linear response as well, but could also be com- 
puted via zone-center "frozen phonons". In both cases, 
the C ss i are the coefficients of the second order expansion 
of the total energy at equilibrium, assuming the ordinary 
periodic boundary conditions for solving Schrodingcr 
equation. In fact, it has already been stressed that set- 
ting periodic boundary conditions is equivalent to set the 
macroscopic electric field E equal to zero; the magnetic 
case behaves quite differently in this respect. 

In MEs a second order expansion of the total energy 
in the sublattice displacements u s , performed with the 
usual periodic boundary conditions, does not provide the 
analytic part of the force-constant matrix C' ss ' appearing 
in Eq. (25). As emphasized above, such analytic part 
implies E = H = 0, while the ordinary periodic boundary 
conditions imply E — B — 0. In fact the real microscopic 
fields, in principle measurable inside the material, are 
E (micro)( r ) and B (micro)( r ). their potent i a l s (scalar and 

vector, respectively) appear in the Schrodinger equation. 
Assuming periodic boundary conditions is tantamount to 
set the macroscopic E and B fields equal to zero. The 
existing linear response codes 3 ' 4 routinely evaluate the 
response to E at B = 0; as for the response to B, the 
needed algorithms appeared very recently. 12 ' 24 . 

Given the above, it is clear that the pair (E, B) is ap- 
parently more fundamental than (E,H). For this and 
other reasons the latter pair has been deemed a "bas- 
tard pair" and an "unholy pair". 8 ' 9 Here instead we have 
good reasons to use (E, H) as independent variables in 
the free energy as well as in the corresponding equations 
of motion, Eq. (20). The key point is that, as emphasized 
e.g. by Landau-Lifshitz, 21 E and H are both longitudi- 
nal, while D and B are both transverse: see Eqs. (2) and 
(15); see also Ref. 22. 

We need therefore to address the following issue. Sup- 
pose we use any electronic structure code which uses E 
and B as control parameters: either setting E = and 
B = 0, or providing the linear response to E and B. 
The quantities entering our main formula, Eq. (25), are 
not these provided directly by the code; we show in the 
following their mutual relationships. 



A. Analytic term 

Suppose that periodic boundary conditions are as- 
sumed to expand the total energy in the zone-center dis- 
placements u s , and suppose that C' ss ' are the second or- 
der coefficients obtained in this way. Which is the rela- 
tionship between these C ss ' and the analytic coefficients 
C ss ,l Setting E = and B = in Eq. (20) at q = we 
get 

fs = ~ C'ss' Us' 
s' 
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(29) 



D 



47T 

IT 



E 



u,. 



(30) 



From the second line of Eq. (30) the actual H field is 



(31) 



hence the sought for relationship between the C sa > and 
the C,,/ is 



47T 



C ss > — C ss > — C s * Moo O ■ 



(32) 



The interpretation of this is pretty clear: the force- 
constant matrix C ss >, computed at zero B, includes the 
restoring forces due to nonvanishing H. This contribu- 
tion must be discounted to get the analytic force-constant 
matrix, which by definition means H (not B) equal to 
zero. 



B. Nonanalytic term 

We rewrite here, for the sake of clarity, Eq. (20) at 
zero sublattice displacements (the q dependence becomes 
irrelevant): 



z*Je 

£ no (%oo 
Moo 



vt 



(33) 



For ordinary dielectrics, the existing codes 3 ' 4 essentially 
evaluate the (Hellmann-Fcynman) forces and the elec- 
tronic macroscopic polarization P linearly induced by E, 
thus providing the tensors Z* and as 



y*] _ df S 

s ~ dE' 



1 + 4tt 



dP 

dE' 



(34) 



This does not apply as such to MEs, because the codes 
implicitly set B = 0, not H = 0. 

More generally, one could envisage a code which pro- 
vides the linear response to E and B: the algorithms have 
been just developed. 12,24 The output quantities would 
be forces / s , macroscopic polarization P and magneti- 
zation M. As for the magnetization, it is comprised 
of a spin (Zeeman) contribution and an orbital contri- 
bution. While the former contribution is in principle 
straightforward, 16 orbital magnetization has been under- 
stood only relatively recently. 25 

We recast the lowest part of Eq. (33) as 



E + 4?rP 
B 



Mo 



E 

B - 4ttM 



(35) 



A linear response code would in principle provide dP / 'dE, 
dPdB, dM/dE, and dM/dB, in terms of which the en- 
tries in the response tensor IZoo could be obtained by 
solving a linear system. 

As said above, the codes existing so far 3,4 only provide 
<9P/<9E at B = 0; whenever the computation addresses 
a ME such response is not simply related to £oo, as in 
Eq. (34). A simple calculations shows that the relation- 
ship is instead 



1 + 4tt 



dP 
dE 



-M 



-v 



(36) 



The correction, being quadratical in a, is definitely very 
small. 

Once the clamped-nuclei response is evaluated as out- 
lined above, we may address the coupling tensors in 
Eq. (33). The field H and the forces are given by 



H = 



Mo 



B - ^alE 



f s = (z;t_ c t M -i a t o)E + c t M -i B . 



(37) 



The last term has a simple meaning: <9/ s /<9B at E = 
is the longitudinal magnetic coupling tensor (analogue 
of the Callen effective charge), not the transverse one 
(analogue of the Born effective charge). The two are 
related by /i^ 1 (by e^ 1 in the electrical analogue). 

The existing codes provide df s /dE at B = 0. In MEs 
this does not coincide with Z*^; the appropriate modifi- 
cation of Eq. (34) is 



_ df s 

s ~ dE 



-v 

o ^OO' 



(38) 



VII. CONCLUSIONS 



Both in ordinary dielectrics and in MEs the zone- 
center dynamical matrix is a nonanalytic function of the 
wavevector q, homogeneous of degree zero in q. We 
have generalized the well established formula for ordi- 
nary dielectrics 1,2,20 to the ME case, where the lattice is 
coupled to both electric and magnetic fields. Our main 
result is Eq. (25), where the coupling appears in a sym- 
metric way. The formula is exact (within the harmonic 
approximation), and is rooted in the formal equivalence 
of electric and magnetic fields in their coupling to the 
lattice in MEs. 

However, the orders of magnitude of electric and mag- 
netic phenomena in condensed matter are not the same. 
ME effects are notoriously small, 5 and the corrections 
to the standard formula for dielectrics — Eqs. (11) and 
(12) — in most cases are expected to be small as well. In 
oxides the electronic dielectric constants are typically 
in the range 2-3, while |^oo — 1| is of the order 10~ 4 , bar 
when close to a ferromagnetic transition. Even the mag- 
netic lattice coupling tensors Q are smaller than their 
electric counterpart Z* by orders of magnitude, e.g. in 
the paradigmatic crystalline ME, Cr 2 03. 7 As for the elec- 
tronic ME coupling, little is known; a pioneering study by 
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Bousquet et al. 16 addressed the Zeeman contribution to 
Ooo, which is of the order 10 in Cr 2 03 (Gaussian units 
are adopted here). More perspicuous effects are expected 
in nonconventional materials, 5 such as those where the 
ME effect can be tuned. 14 

Besides the major result of this work, Eq. (25) we 
have also discussed other issues: (i) The relationship of 
this work to the Lyddane-Sachs- Teller relationship for 
MEs, recently published; 15 (ii) The microscopic origin 
of the coupling of magnetic fields to the lattice, which 
may look counterintuitive; (iii) The relationship to first- 



principle implementations, where in the simplest cases 
E and B (not H) are the control parameters for solving 
Schrodinger equation. 
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